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The single channel autoregressive lattice has been 
successfully applied to problems aa conicene speech analysis 
MismrecOORibion, spectral analysis and noise cancelling. 
More recently the two channel autoregressive (AR) lattice 
has been exploited for autoregressive moving average (ARMA) 
analysis of systems for modeling and identification. 

This dissertation considers the multichannel AR lattice 
when applied to ARMA systems analysis. Constraints on 
Paetee Oaramesers, based on the input output relations of 
the system under test, are developed. The lattice is 
redefined in terms of the frequency domain representation of 
the input data. This proves to be useful because it allows 
the input to be normalized so that the lattice yields a 
consistant set of parameters independent of the test source 
characteristics. Lastly the lattice is redefined in terms 
of correlations of the input signals. This results in a 


computationally andstorage efficient lattice algorithm. 
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In many applications the use of mathematical models 
provides critical insight into dynamic system behavior and 
is essential to many control, simulation and communication 
applications. Models have been applied in the areas of 
economics, neurophysics and astronomy with varying degrees 
of success. Engineers routinely use models when they wish 
to predict or control a systems performance. 

In this dissertation one form of model, the multichannel 
Mietice, is examined in detail in an effort to gain insight 
into its operation and to develop more efficient algorithms 
for its implementation and application to the modeling of 
multichannel linear and nonlinear autoregressive moving 
average (ARMA) systems, based upon a least mean square error 


eeproach. 


fee tik PROBLEM STATEMENT 

The output y(k) of the discrete, linear, time invariant 
system of Figure 1.1 may be described as the weighted 
Summation of the present and past values of the input u(k) 
ema past values of the output y(k). This can be summarized 


Dy the relation 


Pack = 


umd 


3 py (k-n) Cla) 





defining the transfer function 










PACZ)). UZ ) 
Wi 7) eB Cz) 
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ae) Linear Time OR) 


Invariant Digital 
system 






Figure 1.1 General Discrete System 


where an: O0<n <q and one one, uel Pamcdmeters OL the 
system and k denotes the sample index. 

An estimate of the system parameters can be made if the 
output signal, y(k), and possibly the input signal, u(k), of 
a system can be observed. The result is a model of the 
system. There are numerous reasons why this might be done, 
among which are: 

1) The parameters of the system are not known and some 


insight into the system function is desired. 





2) The nominal values of the system parameters are 
known and we wish to detect a change in the system 
Dee rOpmance - 
Equation (1-1) suggests that a system output is predictable 
from a linear combination of past outputs and inputs. 
Because of this, procedures to find these parameters are 
often referred to as linear prediction. 
One way to construct a model is to hypothesize a system 


@f the same form, namely 


y(k) = au(k-n) + b_y(k-n) (1-2) 


en ot ee) 


and adjust the model parameters a. and D_ so that the 
difference between y(k) and y(k) 1s minimized according to some 
criteria. Figure 1.2 shows this method. The symbol "*" is 
used to indicate an estimated value. From Figure (1.2) and 


(1-1) and (1-2) 


Be) = (x) - yk) 


HW 04,0 


Ca a JuCk-n) + 
ale mae) 


gas 


3y¢k-n) - 


Since y(k) is itself dependent on model parameters, and 
a minimum mean square output error criterion results ina 
set of nonlinear equations, the solution of this equation 
-Or model parameters is generally intractable when b_ 


mOmsp is not zero. 
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developed. 


with parameters a. 


parameters, and a synthesis model using a, 






System To Be 
Modeled 


y (k) 


Cie) e(k) 


system Model 


EacuGoute 2 Darect cystem Modeling 


Because of this, other analysis procedures have been 


A very powerful model is the equation error 


formulation, which entails two models; the analysis model, 


' the system. 


The output of the equation error analysis model is 


meacten as 


A q p 
y(k) = ¢f a ju(k-n) vee? 


0 me 


and - which are estimates of the system 


and DD to emulate 


b_y(k-n) Cie) 





where the estimated output is now a function of the system 
past and present input and system past outputs (instead of 
the model outputs). This form is shown as Figure 1.3. With 
the equation error approach a MMSE solution for the model 
parameters is a set of simultaneous linear equations with a 
unique minimum. The synthesis model for the equation error 
approach is of the form of (1-2), using the estimated system 
parameters found from the analysis model. The equation 


Perer Of Figure 1.3 is given by 


R q P 
e(k) = y(k) - yk) = x (a _-a,Ju(k-n) x ° (8 -b_)dy(k-n) 
m=O haw 
(71-4) 
Which is linear in the model parameters. When oe a and 
Bo=b. > e(k)=0. In terms of the z-transform,(1-4) can be 
written as 
fez) = -ACZ)UCZ) + BC z)Y(z) 
SE ee = a) UCz) Gee. 
Bz 
where 
q p = 
Moz = F az and Biz) =1- £ bz™, 
n=0 ei 4 


ale 





u(k) System To Be y (k) 


Modeled 





Pale (S22 
H)AD)yHD 


vk) 


Figure 1.3 Equation Error Analysis Model 


Three cases are examined in this dissertation. 
1) The all pole, or autoregressive (AR), model 
ee = U1 <<q 
2) The all zero, or moving average (MA), model 
BL. = 9, 1l<n<p 
3) The pole-zero, or autoregressive moving average 


(ARMA), model 


is 





an Peeemeceuescome n, 1<n<q 


DL #705 for some n, i<n<p 
The objective will be to develop procedures and algo- 
rithms that allow more efficient calculation of these models 


for application to system performance evaluation and 


ma@fentification. 


Bee OVERVIEW OF THE DISSERTATION 

The following chapters are first a review of some 
contemporary thoughts on linear systems modeling using 
autoregressive, moving average and autoregressive moving 
average modeling, with their lattice formulation, followed 
by some new lattice configurations centered upon the ARMA 
lattice configurations. 

Yule [19] first developed the autoregressive model ina 
paper on sunspot analysis. Because of its wide ranging 
applications, autoregressive modeling has received extensive 
examination since. Chapter II.A begins by reviewing the 
relevant theory of the least mean square solution of the AR 

equations. This solution results ina set of simultaneous 
equations that are the equivalent of the Yule-Walker [9] 
equations. Their resolution requires the solution of the 
feeetx Equation Rb = r, where R is a symmetric and Toeplitz 
correlation matrix. A number of methods have been proposed 
for the solution of these equations. Matrix inversion 
methods can be used, but their computational complexity is of 


concern for higher order problems. Adaptive schemes [18] 
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have been developed that update an estimate of the solution 
as each data point is taken into account. Levinson [6] 
meesented an elegant and efficient recursive method of 
performing this inversion by finding an (nti)-th order 
solution as an adjustment of the n-th order solution. Durbin 
[2] later recognized that the elements of the vector r were 
meso Clements of the matrix R. This was used to improve 

the computational efficiency of Levinsons algorithm. This 
recursive solution can be implemented as a lattice structure 
that has been found to offer an accurate and efficient 
solution to the AR modeling problem. More recently Burg 

[1] has shown that this lattice 1S a maximum entropy 
Petition, proving that the lattice solution 1s optimum. 

In Section II.B the moving average model is similarly 
discussed. The result of a least mean square solution is 
the digital form of the Weiner-Hopf [17] equations and also 
requires the inversion of a symmetric and Toeplitz matrix. 
Perry [16] has shown that this can be efficiently solved 
uSing Levinson's recursion and can also be restructured to 
“form a lattice implementation, part of which is the AR model. 

Sections II.A and II.B deal with single channel (single 
Mput, Single output) models. Section II.C discusses 
multiple input, multiple output generalizations of the AR 
and MA equations and structures that Robinson [16] has 
shown can be solved as a multichannel generalization of the 


Single channel structures. 


ic 





Chapter III is a discussion of ARMA modeling. Unlike 
the AR and MA models, the ARMA model represents both poles 
Ema zeros and so will offer a lower order solution than 
either the AR or MA models for zero-pole systems. In III.A 
the least mean square solution of the ARMA problem is 
developed. A direct application of the ARMA equations 
results in a set of nonlinear, multimodal equations. This 
problem was greatly simplified when the equation error 
formulation [5, 12] was found to result in linear equations. 
Then in Section III.B it is shown how Perry [16] recon- 
figured, with appropriate assumptions, the ARMA problem 
Pemciat it could be solved using a structure identical to 
the two channel AR lattice. This structure offers the 
advantages of a recursive in order solution and results 
in more accurate solutions than batch methods that use 
more direct matrix inversions. 

Lattices have, in the past, used data in the form of a 
discrete time series as a block of data, as discussed in 

Chapter II and III.A and III.B. Griffiths [3] and later 
“Mort [13] developed adaptive algorithms that still use 
discrete time series data but update an estimate of the 
lattice parameters as each data point enters the lattice 
structure. In Section III.C a new formulation for deter- 
mining the lattice parameters is presented. This approach, 
although fundamental, has generally not been exploited, 


Beericularly tor the multichannel case. In this structure 
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the data presented to the lattice is the frequency domain 
Spectrum of the signal to be analysed. The power spectral 
desnity of the error spectrum 1S minimized and the resultant 
lattice is equivalent to that used with time domain data. 
This new lattice structure is useful because it allows the 
spectram content of the input signals to be observed and 
presents the possibility of adjusting the input spectra 
to compensate for a nonwhite driving source. This configuration 
also permits the lattice to be used for filter synthesis. 
Section III.D discusses the lattice ARMA model when applied 
to system identification. The ARMA lattice identifier can be 
used for system performance evaluation. In this application 
the model parameters obtained as a result of a test are com- 
pared with model parameters obtained from a system functioning 
normally and functioning with known faults. This comparison 
can reveal if the system under test 1s operating satisfactorily 
and possible identify a class of faults within the system. 
The lattice parameters fully identify, ina least mean square 
sense, both the input and output of the system to be modeled. 
iron modeling it is generally desirable that the input signal 
be white (have a flat spectrum) so that the signal energy is 
equally distributed across the spectrum of interest. Some new 
relations among the lattice parameters are developed, which, 
with some previously known relationships, reduce the number 


of parameters needed by one half when the input is properly 
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conditioned. Some new methods are then derived that allow 
preemphasis or normalization of the system input to correct 
for a nonwhite input signal. This is important because 
lattice parameters, while modeling the system, are still a 
mm@etion Of the input signals second order statistics. This 
is useful even when the test source 1s purportedly white, 
Sance the lattice must function with estimates of signal 
parameters obtained from finite length data signals. 

Section III.F introduces another application of the new 
frequency domain lattice. The lattice can be used to experi- 
mentally synthesize a filter to meet a desired specification 
by providing the desired response as the input of the 
lattice. The model generated is a minimum mean square error 
approximation of the desired specification. 

In Chapter IV and Appendix A a new approach to both the 
Single channel and two channel lattice is developed. This 
method recognizes that the error signals of the lattice are 
only linear combinations of the input signals and delayed 
versions of the input signals. From this observation the 
equations for the lattice parameters are restructured so that 
ey linear combinations of the input correlation functions 
are needed to solve the lattice. This 1s more computational 
ema storage efficient than other lattice methods since it is 
not now necessary to generate updated error signals at each 


Stage. It retains the advantages of the lattice of orthogonality 


Ike 





and maximum entropy. It offers the new advantage that data 
need not be artificaily windowed to force the correlation 
matrix to be Toeplitz, as is required by Levinson's 
algorithm. This may lead, with some types of data, to a 
more accurate solution. This method is extendable to higher 
dimensioned lattices as are used for nonlinear and multiple 
meput multiple output system modeling. 

Chapter V presents a summary of the new results in this 
dissertation. Included is a discussion of the application 
and implementation of these new concepts, along with questions 


Tor future research. 
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Celio coli MODELING AND LDENTIFICATION 


The autoregressive, moving average, and autoregressive 
moving average models have all been successfully used for 
system modeling and analysis. The AR model is the model of 
choice when only system output 1S available. The MA model 
can give lower MSE if an input signal 1s available, but 
cannot generate an infinite impulse response. These models 
will be examined in greater detail and computationally 
efficient algorithms for their solution will be developed. 
Both single input and single output and multiple channel 


implementations are discussed. 


wee LHe AUTOREGRESSIVE MODEL 

For the autoregressive model the assumption is made that 
the present output of a system can be estimated from a 
weighted summation of the past values of the output. This can 


be expressed as 


y(k) = boy (ken) G7) 


ti mtd 
f- 


@enan matrix form 


y(k) = bi y(k) (222) 
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yi (k) = Ly(k-1) y(k-2) ... y(k-p)] 


If a signal y(k) is to be modeled the prediction error 
can be written as the difference in the signal and its 


estimate or 





e(k) = y(k) = y(k) = y(k) - bly(x) Gan) 
y (k) 7 
271 66 71 
b, dD, 
y (x) 


Pigure 2.1 An AR Prediction Error Model 


Figure 2.1 illustrates this type of system. 


ral 





The mean squared prediction error as a function of the 


weights DL is found to be 


(2-4) 


where in general ale 


R = 
WI 


eGo vecenii > F ap Seco w (ken) 5 
e (x(k) w’ (k)} and e{:} designates expectation. 


In this case, 


0 R (=ptl 

gy 62) gy (TPT) 

R= 

aa Roepe) R (0) 
yy Po yy 

meee (RC) R(p)] 

ey yy yy 


The solution of (2-4) for the optimum weights 1s given by 


(2-5) 


which are commonly called the normal equations. 


foechas aS 
Meimescituted into (2-4) then the minimum error is found to be 


= _ af 
“min 7 yo a Dopt ryy (2-6) 


(2-3) may be written in the z domain as 





=e en) (2-7) 


which is the relationship of error signal output to the 
[moeut Signal. 
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Consider the all pole filter of Figure 2.2 with the 
memos signal u(k) and output y(k). The relationship 


between input and output can be written in the z domain 


a 
a 0 
Ae = H(z) - = 0 - B(z) ee 8a) 
ae biz 
i= IL 
Geen the time domain 
p 
mck) = ajuk) = UD bLy(k-n) (2-8b) 
n=L 





Proms 252) All. Pole Filter 


If the output y(k) of the all pole filter is applied to the 
[emi Of the prediction error model of Figure 2.1 the output 
Seor will be e(k) = ajuk). Miemeweanorenm Function of this 
cascade arrangement from the input of the all pole filter to 
the error output of the predictor will be an: This shows that 


if u(k) is white noise the effect of the prediction error 
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model is, except for a gain term ay to reverse the action 
Of the all pole filter or to whiten it's output. For this 
meson the prediction error model is often referred to as a 
whitening filter or an inverse filter, and makes it useful 
Peeoespecerum analyzer. it should be noted that ay establishes 
a lower bound on the MSE of the AR predictor. If u(k) is a 
Unit impulse, u(k) =46(k), the model may replicate the 
response y(k) perfectly, but can never predict the arrival 
of the input ao (kK) Since the predictor is causal, and there 
is always this minimal error. 

The relationship between an all pole filter of Figure 
2.2 and the analysis model of (2-1) suggests that the all 
pole filter can be used as a synthesis model for generating 


the signal y(k) using ajulk) as the input signal with the 


coefficients obtained from the prediction error analysis 


model. It is only necessary to find a suitable gain term 
ay: Equation (2-3) can be rewritten 
oe . T 
peek) = e(k) + D y(k) (2-9) 


Pemparing this with (2-8b) it is seen that a,julk) = e(k). 
By requiring that the total energy of the output from the 
linear predictor be equal to the total energy of the 


meee tO the all pole filter, a. is specified by 


0 


Z 
ee. ete (k)} 
UU 
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This leads to a z domain transfer function representation 


for the synthesis model of 


ei 
H(z) = 505 Qa) 


The impulse response of (2-11) can be of infinite dura- 
tion. This means that a relatively low order AR model may 
emulate an infinite impulse response system to an 
meeeptable accuracy. 

In the problem as stated so far it has been tacitly 
assumed that the desired order for the model was somehow 
known. In practice this is often not the case. Instead, 
(2-5) is solved for some order and the MSE determined. If 
this MSE is unacceptable, then (2-5) is again solved for a 
higher order. Solution of the simultaneous equations of 
(2-5) requires, by Gaussian elimination, ee + O(p*) 
operations, so that a repetitive solution of the normal 
equations can become computationally ponderous. There is, 
however, a high degree a symmetry in these equations and 
. the exploitation of this symmetry can give a more efficient 
solution of (2-5). 

im A Recursive In Order Solution 

The normal equations (2-5) for the optimum weights 


Gan be written in matrix form 


ZS 





eee cos)! |b R (1) 
ee yy? 


1 yy 
b R.(2) 
Roy () Ry 60) > a 
= (2-12) 
= sis | | b R_ (p) 
Ry (P-L) Ry (0) } , | RO 


The autocorrelation matrix is Toeplitz (the elements along 
any diagonal are equal) if the signals are stationary which 
is a basic assumption of this work. 

Calculation of the autocorrelation functions now 
deserves some comment. Computation of the exact discrete 
correlation requires a summation of infinite duration 


Signals. 


co 


R NOy= 2 k k+N ) 
ol ) ae y(k) y(k+i 


The result is an even function, or 


R_ o(N) = R_ (-N) 
yy NE 


Seema: the autocorrelation matrix 15 symmetric and 
Moeplitz. In the absence of an infinite duration signal, 
the more common case, an estimate of the correlation func- 
tion must be made. If a window function is used and the 
data is taken to be zero outside the window then the 
estimated correlation matrix will be symmetric when the data 


is averaged over the width of the window. If the data is 
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not taken to be zero outside the window and averaging is 
done over a finite number of points, the estimated 
correlation matrix is in general not symmetric unless a 
large number of data points are taken because of end 
effects. In the derivations which follow symmetric Toeplitz 
matrices are assumed. It is noted that the final algorithms 
can be used without assuming symmetric Toeplitz correlation 
matrices, and apparently work just as well. In fact it is 
noted that Griffith's time adaptive technique [3] is 
nonstationary and yet appears to work very well. This has 
not been explained in the literature. If the observed signal 
1s ergodic an estimation of the correlation function can be 
made using a finite interval of data. 

When the data 1s windowed so that the autocorrelation 
matrix is Toeplitz and symmetric, the method used to solve thi 
resultant normal equations is often referred to in the 
literature [9] as an autocorrelation method of solution. 

When the data 1s not so windowed and the matrix is not 
Toeplitz, the problem solution is said to be a covariance 
method. Either method results in an estimate of the true 
values of the correlation matrix, and it is generally 
unclear that one method will be more accurate than the other. 

The Levinson algorithm relates the (n+l1)-th order 
Pemmcion Of a set of equations of the form of (2-12) to the 


n-th order solution. The assumption is made that Sie 


oe 





where the superscript denotes the order of the solution, can 


be found from Sue PD ecodtmuec correction to each term and 


eme new term, or 





pom? | -(n) rg 
ee = |----| + | ------- where ae = i C2= 173) 
0 ent) 5 (n) 
ie _ en 
i tind be parameters en and eee must be computed. 
Permuted versions of the vectors pe and may be 
defined by reversing the order of their elements. 
» on) R sad 
cay we oy, 
__ (a, : n 7 n 
= ° = Ske: = 
2 fyy ; =yy Samet? 
(n) 
Ree 1D 
ig! yy 


The symmetry of the autocorrelation matrix allows (2-12) with 


Capen) (n) 
r 


p=n namely Sos b = , to be rewritten using (2-14) as 


mY 
en) Cn) (n) 
R b = 2-195 
Sa meek 
: ste : : ; : 
The matrix Ron Pe 1s recognized as the matrix Re with one 
~new row and column and may be written as 
eae (a 
(n) : Ee | (n) 
RS JR (=1) Saale 
ont 1) ms SS So eer =o ie 
yy i. ©. | t a - : eng 
Cay) he, C2) RR 2 CO) eee 
Nas | Y Pyy | yy, 
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since the data has been properly windowed to ensure that the 
autocorrelation function estimate is an even function of 
time. The foregoing and (2-13) may be substituted into an 


(ntl)-th order version of (2-12) to give 





Multiplying the matrices and substituting (2-15) yields 


oN = sp Ry Se Cor aisee 
where 
ce Cn) 
(nt1) Rumer) ie) 19) 
k xy a2 (C=ileD) 
Sey a =o 
Mmeetituting (2-16a) into (2-13) yields 
Pe am 
> hs. i fntl) a (2-17) 
0 i 
ee 
iWon eimakie rel Bee From ow it 1s only necessary to 
determine intl) which can be found from (2-16b) which uses 
Royintl) and other terms from the previous solution. It is 


noted by comparison with (2-6) that the denominator of (2-17) 


ZS 





is the MMSE of the (n)-th order solution for the normal 
Pr@ations. Equations (2-13), (2-14), (2-16) and (2-17) can 
be applied recursively, starting with order 1, until a model 
results with the desired accuracy. 

eee AR battice Structure 


(Cn) 


Biemveator ) EnineiemiGiterem Im) che tEransiorm 


domain 


eo, 


ead the reversed in order version b can be written 


362) (2) = 2 7Mgln) (sly (2-19) 


Memsnould be recalled that B(z) is the transfer function of 
the prediction error model. The overbar here is used to 
indicate the reversal in order, and will be used throughout 
Moms dissertation to signify a backward signal, or a function 
- Or martix associated with a backward signal. 

A recursive solution for the transfer function B(z) 
@ameDe found by substituting (2-17) into an (ntl) order 


version of (2-18) 


gintl).y a = reer oS TL) 


Su 





b =-b 
2 3) = ee — 20 -n-1ly cud ae +  ontl) Lore 
0 a 
on) 
: 36M) (5) _ zi, (ntl omnypononel | fa eae 
ip 


and substituting (2-19) into the second term of the foregoing. 
The result is 


Mee) = Be ce) = 2k Bt)? B62) (2) (2+20) 


go) (2) can also be found recursively by substituting 


Mey into (2-19) rewritten for order (ntl) 


Bintl 7) = g7igh™M (2) - ~MtD pl) (2) (2-21) 


Cap 


Sianee B (z) 1s the transfer function of the n-th 
order prediction error model, when it 1s driven by the 
-Signal Y(z) the result is the n-th order error signal, or 


(Cn) 


mez) = BS™ (2)¥(2) (2-22) 


The forward prediction error is the difference in the signal 


and a predicted signal derived from a weighted summation of 


the past n samples of the signal. The backward prediction 


Syn 





error is the difference in the Signal at time (k-=n) and 
its predicted value based on a weighted summation of the n 
Samples to follow. Figure 2.3 isan illustration of elgiglce 


forward and backward prediction. 


y(KkT ) 
vig 
. -1) 
y (ken) ie } e(k) 
e(k) | co) 
y(k- 
k-n kK-entl k-1 k Keen aa 
” 
Forward Prediction 
nN Samples 
Backward 
Prediction 


mreure 263 illtstration of Forward and Backward 
Prediction Using n Points 


By multiplying (2-20) and (2-21) Bye 2) (Zz) ane 
Bi@2) May also be found me Gur save 7. 


gintl)’., = Ef) (5) - x (ntl) -1 E‘) (5) (2-23) 


= \ = = 
Meee (a) = 2) EF) 5) RAF a(n). (2-23a) 
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In the time domain these two equations become 


eee = te Gore e276?) (k-1) (2-24) 
ee = es ety = KBE (xk) (2-25) 
Since the zero order prediction of a signal is y(k) = 0 (no 


prediction at all) the zero order error is 


Bao = 6 ck) = VO (2-26) 


and equations (2-24), (2-25) and (2-26) can be implemented as 


a lattice structure as is shown in Figure 2.4. 


ele 


aS ny Bae 





Figure 2.4 Second Order AR Lattice Model 
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The lattice parameter eee? can be found from (2-16b) 


(as was done for the Levinson algorithm) or can be found by 
Minimizing the expected value of the squared error of (2-24) 
or (2-25). Minimizing square of the the expected value of 


(2-24) with respect to k yields 


ie Gore @o Ge 1) 
(2-27) 


pont) 


3 Cn) 


ie ics =) Cea) 


and similarly for (2-25) 


Pe ree (ics) 


moat) a (le (9293) 
(n) (n) 
ete (k)e Ck} 


These alternative solutions are referred to as the 
forward and backward methods and can be shown [8] to be 
Pemavalent to (2-16b). From (2-20) and (2-21) it is seen 
that the forward and backward prediction transfer functions 
' are equal so (2-27) and (2-28) must be equivalent. However, 
when implemented, the expectations will be estimated from 
time averages and will generally not be equal. This has 
Mean tO two alternative solutions, the first a harmonic 


mean of (2-20) and (2-21) due to Burg [1] 


34 





Me Goa C1} 
(2229) 


D 2 
ze te(™? (one Cem) 
and another from Itakaura and Saito [4] which is the geometric 


Mean of (2-27) and (2-28) 


Meee ie Ca’ Cee1)} 
« int. a C220) 


ies . - . 
ee ones 1) 


Stability of the result is ensured if all the poles of the 
transform domain transfer function lie within the unit 
circle. Markel and Gray [11] have shown that the necessary 


Met t1cient condition for this is 
he? <a C25 


Since (2-30) is of the form of a normalized correlation it 
1s of necessity always bounded by unity, and since the 
magnitude of a geometric mean is an upper bound of the 
magnitude of a harmonic mean, (2-30) and (2-29) are 
guaranteed to result in stable solutions. There is no such 


guarantee with the forward (2-27) and backward (2-28) methods. 





A most important property of the lattice structure 
is the orthogonality of the backward error signals, 


described by 


me Coe)’ (k)} = (2-32) 


Indeed, the lattice structure can be derived [11] starting 
with this property. This orthogonality decouples successive 
stages of the lattice, ensuring that the addition of stages 
in the recursive solution of the filter will not affect 
preceding stages, and permitting the derivation of (2-27) 
and (2-28). Furthermore, [8] orthogonality allows for a 
recursive calculation of the error energy at successive 


stages by 


ee Ge Py) = ete ™ 2) (xy SP ()} = port? 


2 
p(ntl) pf) (5 _y (ntl) ) (2-33) 
Noting that 
Bee = ete) 26 (e)} = RCO) (2-34) 


(2-33) and (2-34) may be used for the denominator of either 


Metop), (2-27) or (2-28). 
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Onee (2-17) is solved, the vector pe may be 


miplemented as the all pole filter of Figure 2.2, or other 
equivalent arrangements, to produce a synthesis model. 

The lattice synthesis structure may be generated by (2=25) and 
rearranging (2-24). 


o ae 2 amr) Cnt LD so 0) 


Coy 4k Ce) C2=)59) 


This equation can be realized by the structure of Figure 2.5 


ae), 


and implements the transfer function 1/B 


CZ 


Ma) —~ e') (x) oD (yy ay(k) 









a2) (K) a aes a) Ge 


mesure 2.5 Second Order AR Lattice Synthesis Structure 


Pee tHe MOVING AVERAGE MODEL 

The moving average model is based on the assumption that 
the present value of a system output can be estimated by a 
weighted summation of the present and past q values of the 


Sieeem input. In equation form this is 


a7 





; q 
y(k) = £ a(n) ulk-n) = a? uk) (2228) 
n=0 
where 
ae 
eet) al) ) ey alg) | 
ultk) = Cu(k) u(k-1) ... ulk-q)J (a7) 


meesnould be noted that the vectors of (2-37) are (qtl) 
elements long. 

Defining the error of the model as the difference in the 
estimated value and the actual value of the system output, 


or 

eek) = y(k) - yk) (2-38) 
the mean squared error E is found to be 

E = ote a- 2 ae r 7 Roe 0) C2= 49) 
where ae is (qtl)x(qtl) and aay is a (qtl) element vector. 
iim this 1s minimized with respect to the coefficients of a 


the result is the normal equations of the moving average 


problem 


2 =e G2 st.) 
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From this the optimum value for a is 


2) = at (2-41) 


(2-42) 






syscem To Be 
Modeled 


Moving Average 


Model y(k) 


Figure 2.6 Moving Average Modeling 


Figure 2.6 shows the block diagram of the implementation 
Seecie MA model. The transfer function of the model is, in 


ihe zZ domain 


-1 
BUC ply (2-43) 


oo 





Mgt e Otten rererred tO as an all zero model. As such its 
impulse response is finite and in all cases the model is 
meapile. 
1. A Recursive In Order Solution 
Pauceeen (2-40)e4s Ssimglar an structure to (2-12). 


R is Toeplitz as is R of (2-12), but the elements of r 
—uu —yy —uy 


~ 


of (2-40) are not also elements of the matrix Rea Sei! 
Levinson's algorithm may be applied. 


Equation (2-40) is written for order n as 


en) Boo) 2 p> 6m) lo sininy 


—uu — —uy 


Grate) ‘Go 


ema the assumption is made that a Pome On Ol aa 


ema a correction vector. 


(2-45) 





(2-46) 


4 Q 





pmaeexpand an (n+l) order version of (2-44) with (2-49) 


poestituced. 
ie (7) (n+1) a” Git) (n) 
aa ore a + y sg 
a ait ts +------- Se eee i atest ea arte 
janie | (ntl) 
5 R (0) E | R(n+1) 
erent | “uu 


When this is multiplied two equations result 


fn). Cn) Gi) end) Gels) eee Cne |) S5- (m) : 
ae = i Sree ib ‘i Pau ~ =uy Oy 
meeey Ge (ntl). itt) Gren 
=) 7] a i erp ie i Ry 69? S = Ryy intl? 
(2-48) 


The first and last terms of (2-47) are equal so both may be 


cancelled leaving, when solved for ae 


-1 
y (ntl) 7 _pin) Gace) ae (2-49) 
- —uu uu 
a. (n)7+ (nn) . 
mending R v is comparable to (2-5), the (ntl)-st order 
—uu uu 
mmmorersression of the input signal u(k). So to derive the 


MA model, the AR model of the input must first be solved. 


ae) 


When this 1s done a new vector can be defined where 


YL 





e(ntl) _ a(n) 5 ntl) 





= = (2-50) 
— —uu =uu 
' 7 : Gin e. 
Mmemmene result of the autocorrelation analysis. Now y Ls 
defined as 
y fay oe (2-51) 


and this can be substituted into (2-48) and solved for go) 


Credle a Gata ee) 
3 ee, = =u . = (2-52) 
abe 9 (ntl) p (ntl) 
uu —uu — 


mee A MA Lattice Structure 

The Levinson solution of (2-49) can be restructured 
Meemeive a lattice configuration. The transfer function of 
the MA model (2-43) can be written recursively in the z 


domain using (2-45) and (2-51) as 


po 1) (n+1)5(nt1) 


C= ie) + g (z) (Pass) 
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where gintl) 


(z) 1s the result of the autoregressive analysis 
Of the input u(k). If both sides of (2-53) are multiplied 
by the input U(z) and then transformed to the time domain the 


result is the estimate of the output signal y(k) for 


“A r “A a = 

ee or yc) & get) (K) (Des) 
where Ge) 1s the backward error signal of the autore- 
gressive analysis of u(k). The desired lattice structure 


MeesnOwn in Figure 2.7. 


PSO) ren of) CK) —(2) 





(k) 


Figure 2.7 A Second Order MA Lattice Structure 
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Recalling that the backward error signals of different 
stages of an AR lattice are orthogonal, it 1s seen that the 


(n) 


MA estimate of y (k) is a linear weighted sum of these 


orthogonal signals. 


See MULTICHANNEL MODELING 

The models previously considered have been of single 
input, single output systems. Multiple input, multiple 
output systems can be modeled as vector generalizations of 
the models heretofore developed. 

1. A Multichannel Generalization 

A single channel output of a Q-channel autoregressive 

model is the weighted summation of the past outputs of the Q 
channels of the system to be modeled. Similarly, a single 
channel output of a Q-channel moving average model of a 
system is the weighted summation of the past N inputs.of.the 
Q channels to the system to be modeled. Either of these 


statements may be written in equation form 


A Q N 
Peck) = »f eee Cn) ac. (ken) C2255.) 
J j21 n=l 741 = 
Q = number of channels k = time index 
n = time delay index } = output channel being 
estimated 
1 = input channel 


uy 





oN = estimate of output channel j at time k 
x, Ck-n) = input to channel 1 at time (k-n) 


d..(n) = weighting factor for estimating the output signal 
tJ see tannewm at time k due to input x (k-=n) 


As an example, for a model with 3 inputs (9 = 3), 


(2-55) would be expanded as 


A N N N 

x Ck) = y d,,¢n) x, (k-n)+ : d,, (n) x, (k-n)+ ; d.,,(n) x, (k-n) 
n=l ig et igi it 

i N N N 

mem= > 3d tn) x.Ck=n)+ © d.,~.(n) x.(Ck-n)+ © d..(n) x.(k-n) 

2 aaa eZ sh ee LZ z. neu a2 3 

a N N N 

Xk) = - d,3(n) x, (k-n)+ - d,,(n) x, (k-n)+ : d,,(n) %q (k-n) 
n=l igus ae n=l 


Equation (2-55) can also be written in matrix form. 


“w 


x(k) = p x (ke) = (x (i pl (2-56) 
where 

x(k)" = [x,(k) x(k) 66. yO] (QxL) 

Mik) = [xy (k) x5(k) oo. x5,() J (NQx1) 

x3(k) = (x, (k-1) x,(k-2) 2. x, GeeNDI (Nx) 
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—11 -12 ~~*° “IN 


5 = dy d,» =e day (NxNQ) 

For 2o2 don 
Pee ice) dd... 02) a4. a>. (N)] (Nx1) 
—1) ay sey 1J 


finemerror or Suchen CStimate is the veetor difference 


ie the actual output vector and the predicted output vector. 
e(k) = x(k) - x(k) = [I - D!] x(x) 
he prediction error covariance matrix is a (QxQ) matrix 
q T 
m= efe(k) e (k)] (2-57) 


If the trace of the prediction error covariance matrix 1s 


meeimized the result is 


RD-r (2-58) 
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where 





: i 
Be ce Xeec) RMCK)} = € x5 () J 
Recs 
S16 
= (2259) 
Ree. 
=o —C 





M——eMemarrix Of aucocorrelation and cross correlation matrices 


and 


K(k) J = 





(2-60) 


"1S a cross correlation matrix of size (NQxN) and elements 


ae Ci) 
io] 


jx, | 
— 


(2-61) 
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As imposing as this may seem, it 1S no more than a 
vector generalization of (2-12) from the AR problem and 
(2-40) from the MA problem. Each correlation coefficient is 
replaced with a correlation matrix in the multichannel 
generalization. 

It will be useful in the sequel if the two channel 
lattice is examined additionally. Equation (2-56) expanded 


mer two channels is 


T 
0 

rr 

2 
T 





N d,,(n) “am x, (x-n) 


d,,(n) dy y(n) | gy Sheol) (2-62) 


Me two channel error is the vector difference in the signal 
to be estimated and the estimation. 


x, Ck) x, Ck) (x) (x) -x, Cx)) 
e(k) = - - (2-63) 


SSS %, Ck) oS) 
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Equation (2-58), the result of a minimization of the error 


equations, is in two channels. 








doo Px ames et x er 
| al Da, 
where 
R (0) a ein (1-N) 
re »& . 
sg i 
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—X .Xx 
i) 
R Ge) Ry (0) 
alas a] 
R ” i) 
ca ey 
eS cao 
Ry - (CN) 
Cae 5 
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hieectanestoutie Solution of (2-58) is the form of 
the (NQxNQ) correlation matrix R. The matrices on the main 


diagonal are all correlation matrices of the form R, e e and 


so are symmetric and Toeplitz. The off diagonal ai 
matrices are cross correlation matrices that are Toeplitz 
but not symmetric. Symmetrically opposed matrices are, 
however, transposes of one another, with the result that the 
matrix Ris also symmetric. This block Toeplitz form is 
amenable to a block form of Levinson algorithm solution. 


The vector generalization of the AR lattice is shown 


in Figure 2.8 and is described by the equations 


a 
eGo ee ck) = SPP? SS) (e221) (2-68) 
ne) =a) Aigo) Gp 
@ Goes 2 1) = kK oe) (2-66) 
Mmeemcne initial conditions 
eV (K) = 269 Ck) = x(k) (2-67) 


The signal paths of the single channel AR model have been 
replaced with parallel vector signal paths, adders with 


vector adders and multipliers with matrix multipliers. 
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Fagure 2.8 A Two Stage Multichannel AR Lattice 


(n) (n) 


Eteeice predictor coefficients K and K may be 
also found from vector generalizations of the single channel 


case using 


ait T 
at} 2 ae re (2-68) 
ent) ) Ca) tn) 
K as A (2-69) 
ea) Ge Gaye 
Mee ete (k) 6° (k-1)} (2-70) 
a) aD Gee ceap CRO 
Mee = ete "7 (k) ec’ (k)} = OP a Seger enuian 
7 
=) 2G) eee Sea eG 
Mee = cite (k) a’ (k)} = OP tc = ear 
(2-72) 


OL 





As in the single channel case, the optimum predictor 
for all orders less than N are found with the N-th order 
solution. Also, orthogonality applies to the backward 
channel error signals, or 


: oa 
efe-(k) @) (k)} = 


and successive stages are thereby decoupled. 
For a multichannel MA model equations (2-65) to (2-72) 


are used with a vector generalization of (2-54). 


= nN fh 
me) (x) 2 eG.) i: avo) ah Re, 


(k) (2-73) 
The error of the MA model is expressed by 

e§™) (k) = y(k) = y (kK) 

—O — — 


meorcan be minimized with the result 


fees s(n) ele" (k) y' ¢k)) (2-74) 


ro) 
tI 
lrul 


Figure 2.9 illustrates the resultant structure. 
In the next chapter the multichannel generalization 


is applied to ARMA modeling. 
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u(k) 


Pasure 2.9 


A Two Stage Multichannel MA Lattice 
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Mew eee UTOREGRESoLVE MOVING AVERAGE MODEL 


Pie ibe eeREQUENCY DOMAIN 


While the MA and AR models are suitable for many appli- 
cations, modeling some systems can require a very high 
order model to obtain suitable results. Consider the 


geometric series expansion. 


iL 


Cie az t 


This shows that a single pole is equal to an infinite number 
of zeros, and conversly a single zero iS equal to an 
infinite number of poles. Thus a system with a single pole 
will require an AR model of infinite order, as would a 
system with a single zero require an infinite order MA 
model, to be precisely represented. The ARMA model emulates 
both poles and zeros, and so can model pole-~zero systems 
without resorting to infinite order solutions. 

This chapter reviews contemporary ARMA lattice struc- 
meees, then develops some new formulations using frequency 
domain data and considers their implementation as a 


performance evaluation tool. 
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A. THE ARMA EQUATIONS 

Hommodel sa system with both poles and zeros it is 
Macuitively appealing to construct a model that also has 
both poles and zeros. As was earlier discussed, the equation 
error model will result in linear equations as the MMSE 
solution. If the condition that the order of poles equal 
zeros is imposed the equation error model output estimate is 


described by 


N N 
5 by (k-n) (3-1) 


y(k) = ¢f a u(k-n) ~ 
n=0 n 


ty7 C0 )a7G01 =| 
2 | 

where vectors y(k), u(k), b and a are defined and dimensioned 

as in the AR and MA discussions of Sections II.A and II.B 

and y(k) is the ARMA output estimate. This constraint is 

not limiting since some of the coefficients of the numerator 

.or demoninator can be zero. Minimizing the mean squared 


error with respect to the parameter vectors a and b yields 


SS 





| . Be aaa 
(0) .6. RY C1-N) Roy CL) Ry yO) ++. RG] fb | PRY 




















” i C8. b 
R__(N-1) RY (0) aa ma ) N] = 
Riy6O) see RyyCAN RCO) RCD)... RUC) | fay |R,,(0) 
° Riyy (4a) ere) eee CO) ae: igo de=hie =a 
° ° re ° e F 
: : | ; : 

Ry (NFL). Ry) FRAG RA GH).. RCO) fay) | 
(3-2) 
fen a MMSE of 
E Tat ae le, : 
en = Role - Lb opt > ae yy (3-3) 
rn ND) 
ee 
Ree Go) 
uy 
R (N) 
uy 
_ = e | 
Memeacion (3-2) can be readily solved for opt and A opt DY 


matrix manipulation but again a recursive in order solution 


mean be derived. 


fees RECURSIVE IN ORDER ARMA SOLUTION 
Mine (C2Nt1)xC2N+1)) correlation matrix of (3-2) is 
neither Toeplitz nor symmetric, but is nearly so. The 


offending elements are those associated with the ay term. 
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If it is assumed that a, can be found by other means, (3-2) 
can be restructured so that the correlation matrix is block 
Meenlitz by subtracting the a, terms from both sides of the 


matrix equations (3-2). 








yy? a LwlR, C0 «8, (2-90] [b,] By) 
R, (WeL) R60) La gONED) R60) a C0 
Ray? + -Fay(DW Ry) «Rag - Ray) 
Ryy(NHL) Ry (0) yA yO Ruy Rg 
(3=4) 


Since ay is just the DC gain, and can be easily calculated, 
this is a reasonable assumption. Equation (3-4) can be 
rewritten using the notation of (2-64), as was used in the 


memeewssion of the two channel AR lattice in Section II.C 








R » | r ie i‘ ri 
al | | —yy ya 
Ser (3-5) 
R 2 iG i -a,| 
gop ee 
where the vector a contains the N elements Ay ay Of che 
vector a. Here use has been made of the identity Ray =Roy: 
This is compared with (2-64) with x, = y and x, = u 


3k 2 


Sa! 














To solve (3-5) using the algorithms for the two channel 


lattice it is only necessary to make the substitution 


(3-6) 





after solving for the vectors es iy (CNewle1G tele) gaedyel ley rele\slauele 
Since ay is the DC gain of the system to be modeled it can be 


Bound using 


Mey Ck ck) } 
“5 = BIGGS RUG a ar.) 


ime An ARMA Lattice Structure 
Since the equation error ARMA model can be solved 
in terms of a two channel AR recursive formulation, it is 
Meeediiticult to modify the AR two channel lattice to 
Pepresent the ARMA system. 
The ARMA error can be written as the difference 


between the output y(k) and the estimated output y(k) or 


ate! 





=) b 
@, (Kk) = yk) Ly" Gx) a x) LE (3-8) 
from the error equation model of Figure 1.3 and (3-1). The 
Peeerror can be written by substituting (2-62) into (2-63), 
meen y(k) for x4 Cx) and u(k) for x Ck) ; and transposing, 


so that 


A. 


Ak 

y(k) | deed 

e(k)? = wa : Cy) ah) 7s 212] (3-9) 
—21 —22{ 


fa (i) aGouk 
VY u 


where u'(k) is an N element vector containing the elements 
feet )...u(k-n) of the vector u(k). 


hae Jets post Muetiptied by a vector wy where 


iL 
p = - Ceo) 


jihe result is 


T a 
e (k)v Le (k) e (k)] 


Ie 
= y(k)-apu(k) - Cy(k) ut (k)] my Cares 
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file is exactly the ARMA error of (3-8), Furthermore, the 


mean squared ARMA error is 


i 


{es (k)} ye e {e(k) e(k) by 


eek) <2. (Kk) e (k) e (k) 1 
y yy Ns u 


| 


[1 - a,J Ee 
amtk) er-(k) ek) eck) =a 
u y u u 


= bea Cap 


em@en 25 the forward prediction error covariance matrix of 
the two channel AR model. If (3-12) is minimized with 


respect to ay the result is 


efe (k)e (k)} 
ay = — Cae) 
ele (k)} 


Meech 1S an alternative solution for ay 

Pomeonstruce an ARMA lattice it 1s only necessary to 
derive the ARMA error using (3-11) from the AR lattice. 
A second order lattice model is shown as Figure 3.1. 


From (2-65), forward prediction in the two channel AR 


lattice can be written 
2.0) (n+1) 2 Cae) | ™ 1c te (n+1) a dc1) | 


2, (8) : 2,00] a ae E,Ge-1)| 
7 C3 ie) 
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From this equation 0 (K) and a 


(n) 2 nd) (n+1)— zfn) Cee (n) 
=o (k) = ey (k) + Kaa eos (k- -l)tk,y oe (k-1) 
(3-15) 
(n) ox) - _(Atli iy) fl oa S SOs L+kset Ss 87) (K- 1) 
wt u LZ 
(3-156) 


Boammensn(o=15) and (3-16) describe the synthesis 
structure that can be used to implement the ARMA model as 
S@own by Figure 3.2. For an N-th order realization an impue 


(N) 


meeneeded for a Ge). Proms (3—=11) 


= 2, (x) + 4 ai ld Ges (3-17) 


BON) 
Cu 


and if the ARMA model is an accurate representation e, fk) 


Will be small so that (3-17) can be rewritten 
Con =na en (k) (3-18) 


[Mies provides the input to Bp ShoOwmein Figure 3.2 


Meet Tt tCE ALGORITHMS IN THE FREQUENCY DOMAIN 
mr the foregoing it has been tacitly assumed that the 
data available for analysis is in the form of a finite 


length time domain signal. This may not be the case. The 
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(k) can be obtained. 
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data could be in the form of the Fourier transform of the 
time domain signal, or it may prove advantageous to 
process the data in that form. Makhoul [9] has shown that 
linear prediction can be done equivalently in either time 
and frequency domain, but there has been in the literature 
no exploitation of this for the lattice algorithms. Thus 
it seems useful to examine the lattice structures when the 
data 1s represented as a frequency domain spectrum. 

eel ne ingle Channel AR Lattice 

The time domain lattice was formulated using the 

time domain equations (2-24) and (2-25). With frequency 
domain data available the single channel AR lattice can be 
Formulated using (2-23) and (2-23a). Its structure is shown 
as Figure 3.3. A unit time delay becomes a multiplication 
=] W 


by e with the input signal represented by the continuous 


spectrum Y(w). Note that in the backward channel the signal 
a Cu) is a delayed version of Oe Gane The frequency 


domain error of a zero order prediction is 


(0) (3) 


E Mii sale! Cin) =) ecu) ‘Gers 


From the forward error channel recursion (2-23) the magnitude 


seuared error ef order (n+l) can be written as 
et ee Ge ae VE ™ Coyare’™ cw) 


(Atle) (yy] (3-20) 
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Peeure 3.3 A Second Order AR Patan ree in ener requency 


where E*% indicates complex conjugate. The total magnitude 


squared error can be found as jhe ineegral of (€3=20) 


oO 


2 
Mee (oylaa = ¢ JE" (a) | ‘a 
moe ee Gane Cuoidaty TES? (aE (w) del 
Cer aC). 2 
ik pO |B ca) rol (3-21) 
Since XZ" = [X#Z]* the bracketed term of the right side of 


(3-21) can be rewritten as 


695 


Je 


Pf CES ™ (mE (ado + ¢ CES®? cw E(w) Jaw] 


af 
oO 


Goatees: ts Ge (ad) de 


OO 


=(n) 


(3-22) 


Since Y(w) 1s the spectrum of a stable and causal sequence 
Siem Re{Y(w)} is an even function and Im{Y(w)} is an odd 
Function. This is also true for the error signals BY Gap 
and a). The result of the integrations of (3-22), must 


therefore be real and the two terms of (3-22) are equal, or 


[ee] 


Cn) yp)? aps cone ee ref) Cw) ES) (Hy) Jaw 


=e CO =e OO 


ete Ge Cajnide = 27 Ce cE (a) Jadu 


— CO 


C3=23)) 


Soweto 5-25) for the middle term of (3--21) 
and minimizing the result with respect to rahe yields 


ff res™ (yy EM) Cw) Jdw 


ntl ey for a continuous 
«S ) = —---——___———— Ss spectrum Yu) (3-24) 


f |e ean cia 
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or equivalently 


f, 


T 3 
f te Ga (w)Jdw for a periodic 
Mmontl) _ -7t spectrum Y(w) _ 
> : with period 27 S208) 





a 
@ fee) in 


or 
mel ( zs . 
5 pit) (k_ye() (K_) for a discrete 
m1) =e (oul ee mT periodic m-point 
ik = a spectrum yc Ky (3-26) 
eG “m with period mT 


If the discrete correlation theorem and Parseval's 


theorem are applied to (3-26) the result is 


e(nt1) 2 fe age (k-1)} 


: Se ty, (3-27) 
efei™ (k)} 


F 


Mummeh is equivalent to the forward solution (2-27). 
The same procedure can be applied to the backward 


Pmemer recursion (2-23a) with the result 


f ce 2? Gone (w) dw. 
eer) es Seva COnertnuous 
k = spectrum Y(w) C320 


Pas? any ire 


= feo 
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1 Gy = 
ee - ne (w)E (w) Jdw Meese icatic 
k = spectrum Y(w) C270) 
p fe ™ (yy | Pau with period 27 
-T 
m=-1 = 
z ge Cx) Tor a discrete 
Gas 1) k=0 ny th periodic m- 
k = C3 =310)) 
m=. ey 9 point spectrum 
EE Ga? y(S) with 
k=0 mT 


period mT 


which can be shown to be equivalent to (2-28), the backward 


equation. 


—(n) 


eae” (R21) 


(ntl) z , (ntl) 


= ke Coa) 


Z 
e{ei™ (k)} 


2 lhe Two Channel ARMA Lattice 
The multichannel lattice can be similarly developed 


using frequency domain transformations of the time domain 


equations. Doing so, the two channel lattice equations 
become 
| Ee ) 
mes) eet? Gm) = | 262 (3232) 
= a CLG), 
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-] Ab 

1 ae 2 pe ie yeas 
ee) Gay Gas 
K P A (3-38) 

m1) en) Gael) 
oa Ce a (0) =K Is 148) (3-39) 
Ga) ey hoe) Ea) 
EB (w) = E (w)-K E (w) Se) 


where the limits of integration depend upon the periodicity 
of the available data. If the data is in a discrete form 
the integrations are replaced with summations. 

It 1s important to note that when the data is ina 
discrete form the result of a correlation obtained for a 
time shift other than zero will differ when calculated in 
the frequency domain from the result normally obtained from 
time domain data. Consider first two data signals y(xk) and 
Mees DOth defined for 1<k<N, and k an integer. The 


correlation of delay 2 is obtained by the summation 


N- 9. al 
Roy 6%? > a y(k) x (k#&£)| y= 
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and has the effect of windowing the data so that data points 
outside the window are ignored, If these data signals are 
transformed to the discrete frequency domain with a digital 
Fourier transform the resultant spectrum is periodic, which 
implies that in the time domain the data signal has also 
been made to be periodic. A correlation obtained in the 
frequency domain is now 


TT 


Dae ; iw 2 
Riy(®) = Gof Y(w) X(-w) eit dw 


Y coma 


which when transformed back to the time domain is equivalent 
Me, 
2 N 


PE) <a) oF z y(k) x(k+2-N) Vou 
=1 k=N-2£+1 


— 


This is because the effect in the time domain of the fre- 
quency domain time delay elwd feed eCmrenlar Shaitt of the 

data in the time domain. As was discussed in the pre- 
-ceding chapter time domain data would normally use a linear 
shift and then be either windowed or would be considered of 
M@itanite length. Either method introduces it's own dis- 
tortion into the estimation of the correlation. The circular 
Bmitt of the frequency domain correlation is a third and 
Beerterent kind of distortion. This is important when the 


Number of data points is small, because a bias will be 


introduced into the estimate of the correlation. 
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ee oyolteM EDENTIFICATION WITH THE ARMA LATTICE 

The problem of system identification concerns the 
determination of a mathematical model of a system from its 
mete Output relationship. An impetus for this research has 
meen the issue of fault detection within a system. Ifa 
model can be constructed of an operating system, either on 
line or under test, the model parameters might reveal if the 
Beotem is functioning normally and if not, what class of 
component is faulty. While the ARMA lattice has charac- 
teristics which make it particularly attractive for this 
application, particularly its robustness, there are some 
characteristics that still deserve attention. 

In system identification the lattice model would be 


applied as shown in Figure 3.4. The lattice parameters pons 


go and a, that result from a test can be compared with a 
dictionary of parameters obtained under normal and faulty 
Conditions. Using lattice parameters directly, though, 
meets in more model parameters than if the vectors a 
meee were used for identification. Identification using a 
wmaeb uses (2ntl) parameters for an n-th order system while 
lattice parameters result in (8n+l) parameters. Therefore 


it is beneficial to reduce the number of lattice parameters 


that must be examined in the identification problem. 


Ud 
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meeuire 3.4 ARMA Lattice Used for System Identification 


1. Restraints on ARMA Lattice Parameters 
ARMA lattice parameters are a function of input 
Signal second order statistics. The input signals u(k) and 
y(k) to the ARMA lattice are related by the transfer function 
of the system under test. This suggests that the ARMA 
lattice parameters are interrelated. To examine the pro- 


blem the equations for the lattice parameters will be 


jes 





expanded in terms of the system input u(k) and output y(k). 
The vector defining the initial conditions for the ARMA 


Hattice first stage is 


ye) (9) = 
mck) = =e OS (k) Cu) 


u(k) 


The forward and backward prediction error covariance 


matrix become 


T 
2 

vac) Cie UK ie) 

SG (3-42) 
y (k) uk) woe 

ama the matrix ys becomes 
(0) COD ja easlelO 
A = e[e’-‘(k)e* *(k-1)] 

yt) yCk=1) pes) 

= € sats) 


ulk)y(k-1) wO)uCe=D | 
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Now consider the result when the signals u(k) and 


y(k) are from a system as shown in Figure 3.5 






White Notse 
source, 
PSD= ony 









Figure 3.5 System Excited with White Noise 


From linear system analysis theory several statements 
can be made about the second order statistics of y(k) and 


u(k), assuming that u(k) is white noise. 


a for L=0 
eeutk)ut(kt+2)} = Royl® = Coe) 
0 cer 220 
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etutkiyCk+s)} = Rk) = (3-46) 


Bey iiejuckte)} = Roy” = Re 8S 


(3-47) 


h(0) =a (3-48) 


Wieser emel te ecouocan be substituted into (3-42) and 


+3). Thus 
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With these, the lattice parameters can be found 





a 2 ; 
O =o el Ree Ci) ¢ 
el) 2 plo Co" 1 H a yy 
F = ia pent o) |ao2a. RCO! o2ntl) 0 
= ate O yy eee 
ee Z 
S-(R(1)=a,on(i)) 4 
i 10 VV u 
Gs ee Gee Condi )-ak(1)) 0 
= br ayy. O yy 


. a -o-a, | fR,(1) ohh) 
fae =. p80) 00) ~ i. : yy, 
= = pet (P62) |-g2a. RB (0) 10 0 

= moO yy e 
ch) ee 
1 BE Li 
a | ' (3-52) 
Dette ) ote ~G a,h(1) 


mao o-oieeanmad (8=57) Some restraints on the 


lattice parameters are evident. 


we CL) 
fia = Koo = 9 (2260) 
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Ko) > a gka 4 (3-54) 

—({ 1) = —(1) 

Ko 2 = ayky 5 (3-55) 
en. —C 1) =U alee) 

Ki = Ki} + Kyo (3-56) 


Became ens be =seen that 14 1S not necessary to 
calculate all eight lattice parameters. Only three, plus 
the DC gain a,» are needed. Alternatively, the parameter 
matrices can be calculated and either (3-54) or (3-55) used 
for an alternative equation for calculation of a). 

Equation (3-53) can be interpreted in light of the 
Beecice structure. ai and ag are the coefficients that 
eme the predictors of the future value of u(k). Since u(k) 
Mosuncorrelated from sample to sample no prediction can be 
“made and these must be zero. 

Some of these relations can be extended to later 
Merees. The forward error vector of the first stage is 


FP ies aly. . 
* , ty ieee ce) uC] 


e (k) = (3-57) 
uCk) 
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because of (3-53) while the backward error vector is 


P,ty(k-1), y(k-2), u(k-2)} 


e Ge=) ss (260) 
a aKa 2))"s u(k-1), y(k-1)} 


where By indicates a linear function of the arguments. 


mn? 1s now defined as 


abt) 
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a) Gly 
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0 0 
because of the independence of samples of u(k). When this 
Meeapplied to (2-65) for ee the result 1s 
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By induction this argument can be extended to 


arbitrary order with the result 


Gao. .(n) _ a 
Ki9 = Ky =Ormety 2=. 0.5, ley 


when 
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BI u 


u 


Cs G7} 


Using (3-61) the forward error vector can be written 


for order n 


a u(k)+F, fy(k-1),.. 


syy(Kon) UCKeL) y+ suCeon)} 


aCe | 


(3-62) 


mage tnis can be used to find the forward error covariance 


matrix 
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etu(k)aju(k) tu(k)-F, ty(k-1),.. 


ie ee 
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(3-64) 


TE this is now used with (3-59) of order n to find 


the backward parameter matrix the result is 


[2 i) Co (n)] 
= Cu ead At 9 
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Ko =a S Ki5 (3-66) 


For the single channel AR lattice the forward and 


backward coefficients were equal, but for the multichannel 


lattice it is known that the forward and backward parameter 


matrices are not equal. 


determinants 
matrices are 
eeeerminants 


Matrices are 


Nuttal [7] has shown that the 
of the forward and backward error covariance 
equal, and from this it can be shown that the 
of forward and backward lattice parameter 


equal. 
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Tt is knéwn that P * = SOME rollows that 
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pet (p69) = pet (B®?) (3-67) 
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jm is known that 
Sree) = Der(AyDet(3) = Det(B)Det(A) 


Meo that the determinants of (3-68) and (3-69) can be written 
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a 
Minich shows that Dace et”) = et 


) if Det(P 


(1) ) pet (P 


Maeyiew of (3-67) it can be inductively reasoned that 


pet(p'™) = pet (PS?) 


(3-70) 


Now the determinats of the lattice parameters defined by 


=o8)and (2-69) can be written 


-1 
pet(x'n*1)) = pet (P’™ Det (A 
(ntl) (n)7t 
Det (K ) = Det(P )DetCaA 
and since Det(A) = Det (a!) 
pet(K°"*+? 7 ee ee 


A summary of these restraints 1s presented as Table I. 


TABLE I 
SUMMARY 
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Bet (Kk : ee Det (K'™) 
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Zoe Unioueness of the Lattice Parameters 


A particularly troublesome consideration in using 
ARMA lattice parameters for system identification is the 
lack of uniqueness of the lattice parameters for the system 
under test. This is evident because the lattice parameters 
are not just functions of the system under examination, but 
also of the test signal used to drive the system. Thus for 
different test signals such as an ensemble of finite white 
hoise input signals, an ensemble of lattice parameters are 
obtained. Interestingly, experiments have shown that the 
transfer function coefficients calculated using these lattice 
Parameters will generally have less dispersion than the 
lattice parameters from which they are calculated. 

Thus, the selection of a test signal influences both the 
@eeuracy of the model and the expectation of the lattice 
Parameters that result. The mean square value of equation 
error of the ARMA model can be expressed in integral form as 


TT 


Es = weed @ 2) BeCad ei area 
—- TT 


The equation error can be written in the z domain as 


meee = [B(z)H(z) = Az) lUCz) C725) 
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PMpstituting (3-725) into (3-72a) results in 


T - ° ‘ - 
B= af [[BCet Het? )-aced Oh) YA ued?) | Pawt 
—-T 
(3073) 


Equation (3-73) shows that the power spectrum of the driving 
source is a weighting function for the resultant modeling 
error. To obtain a model that is equally accurate throughout 
the spectrum, a test source must be used which has equal 
energy throughout the spectrum. An obvious source with a 
uniform power spectrum is a white noise source. Besides 
giving equal weighting to the error spectrum, white noise 
involves, as has been shown in the preceeding section, 
Simplified calculations and a simplified lattice structure. 
This, however, is not the total answer. When applied to 
system identification, the correlations and cross correlations 
of y(k) and u(k) generally will not be available. Instead 
these quantities will be estimated from finite sample 
.functions of y(k) and u(k). Even though these signals may 
be ergodic, the values obtained are only estimates of the 
correlations. 

Experiments with a white noise source have shown 
that different ensembles of input signals will result in 


different lattice parameters for each member of the ensemble. 
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In order to improve the ensemble averages obtained for the 
lattice parameters, some form of input preprocessing can be 
Bertormed to whiten the input or equivalently flatten it's 
Spectrum. 

Consider first Figure 3.6a, which is a modification 
of the system of Figure 3.4. An adjustable filter with 
transfer function H(w) has been added to the input of the 
system under test. The transfer function H(w) oa be adjusted 
to made U'(w) white if the spectrum of the noise source is not 
white and has not zeros on the w axis. If the system under 
test is linear, then the lattice filter derived from Figure 
3.6b 1S equivalent to the one derived from Figure 3.6c. 
Meeure 3./ 15 a block diagram for calculating the frequency 
domain ARMA lattice model with four preprocessing techniques. 

a) Preprocessing in the time domain on both the input 
and output signals of the system under test. 

b) Preprocessing in the frequency domain on both the 
input and output signals of the system under test. 

Since time frequency domain operations are equivalent 
they can be interchanged. For comparison they are considered 
separately to determine if there are any computational 
advantages of one over the other. It should be noted that 
as long as equivalent operations are performed on both input 
and output, as discussed earlier, the resultant lattice model 
fee Still model the system under test, though the specific 


Battice coefficients may differ. 
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Figure 3.6 Lattice Input Preprocessing 
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Figure 3.7 Frequency Domain ARMA Lattice with Preprocessing 
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A random noise generator can be represented as shown 
moms s.6c. It this generater provides the source signal 
See btgure 3.7, then preprocessing to make U'(w) white is 
equivalent to the multiplication in the frequency domain of 
UCw) by Se), assuming that the inverse exists and is stable. 

In the discrete n-point frequency domain U(w) can be 


merrtecen as a vector 


U J Cai=7 tae) 
Peerereim its discrete form can also be written as a vector 
mee= (GCG. G..... G J (3-74b) 


mietiplying by Seamer can be expressed in the discrete 
Meeduency domain as a multiplication of the vector U by the 


Square matrix H, namely eee where 

me= Eh..] and eee 0 s@ope) a ay (3-75a) 
and 

ime. «= «UC. nue ge lr =eiw le alae ox > ye 1 (3-75b) 


where U; is the expected value of the ith harmonic of the 


input signal U(w) taken over an ensemble of the input signals. 
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Poe eeomecmmactmi SCH) should have no poles or zeros 
on the unit circle. G(w) can be known apriori (or estimated) 
from an autoregressive analysis of the signal generator 
output, or by averaging the spectra of several (possible 
overlapping) segments of the generator output. Equations 
(3-75a) and (3-75b) are useful when dealing with an ensemble 
of input signals. Alternately for a single member of the 


ensemble 


Wwe U- . for i= 1, 2, ..., n (33) 5e) 


can be used. This is equivalent in the time domain to 
exciting the system with a unit impulse. 
paommerwemire 3.6c UCw) = N(w)Gtw). From Figure 3.6a 


the inputs to the lattice filter are given by 
Uw) = Ulw)HCw) = NCw)GCw)H(w) = Nw) 
when 


me) = G -(w). Also, 


fms) = U'(p)TCw) = NCwIT Cw) 


and the transfer function derived by the lattice model is 


os = Tw) 


oe 





From Figure 3.6b U'(w)=U(w)H(w) and Y'(w)=YCw)H(we)=UCw)TCw)H(w). 
ie transiter function for the lattice is 


Gas) 


wey 


which is equivalent to that of Figure 3.6a. In Figure 3.6b 
the spectra presented to the lattice filter are not functions 
of G(w), and the coefficients obtained can be used directly 


to identify the system under test. 


mee EREQUENCY DOMAIN ARMA LATTICE SIMULATIONS 

To verify some of these concepts a FORTRAN program was 
written that implements the ARMA frequency domain lattice 
with preprocessing introduced as shown in Figure 3.8. A 
Hamming window can be selected for preprocessing in the time 
domain. The selectable frequency domain processing imple- 
mented is that of (3-75c). The spectrum Y(w) is divided by 
the spectrum U(w) to give Y'(w), the new system output 
spectrum. The input spectrum U(w) is replaced with a unit 
magnitude real spectrum. 

When preprocessing is not desired then the inputs to 
the ARMA frequency domain lattice are changed to U'(n)= 
mem) and Y'(n)=Y(n). 

For comparison a program was written that solved the 
me normal equations by Gaussian elimination. Listings 


of these programs are available from the author. 
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More than twenty different systems, of up to fourth 
order, were used as the system under test. The number of 
data points used for a single analysis varied from 32 to 
4096. The number of points was always an integral power of 
two because of an FFT routine used in the frequency domain 
Hattice. 

ieetilustrate the ability of the frequency domain 
lattice to accurately model a system two examples, typical 


of all the problems exercised are shown here. System l, 


esis, Ie 
SleS IM de 


ie ticsecm hime tlon CoetTricients 


Numerator Denominator 
Peay) = 0.25 Bolje= =1.14 
meet) = 0.35 BG 2) as Slaeae nice 
meer) = 0.245 Bose = —0.6 49 
Pees) = 0.0 ey ee Es Ort 
ze!) = 0.0 
Zeros Poles 
Real Imaginary Real Imaginary 
-0.70 Cm 70 0.510 Ger50 
-0.70 -0.70 Gyo 8 -0.50 
G07 0.90 
Ope, -0.90 


characterized in Table II, is a fourth order low pass system. 
Figure 3.9a and 3.9b show the magnitude of the transfer 
function of the model obtained by exciting the system with a 


Bangle (256 point) sample of white noise from the Gaussian 
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elements and the frequency domain lattice using no pre- 
processing. As has been observed in the past [15] the 
Tattice results in a closer model of the system, Figure 3.9c 
and 3.9d show the pole zero plots of the resultant transfer 
mimetiLons. 

Figure 3.10a and 3.10b are the transfer function magnitude 
found with the same methods using 4096 data points. Here the 
lattice model is still a more accurate model than the batch 
method. 

While in most cases the lattice offered superior perfor- 
Mance to the Gaussian elimination method, System 2, Table III, 
is a counter example that shows that thisis not always the 
case. In this example the system 1s a second order band pass 
fmeeter. Figure 3.lla and 3.11b are the transfer function 
magnitude plots found using 128 points, and Figure 3.1llc and 
3.l1d show the pole zero plots. Here it is seen that the 
Batch method produces a superior model. Figure 3.12 shows 
the same system analysis using 1024 points, and both methods 
accurately model the system. 

To determine the ability of the input preprocessing to 
compensate for nonwhite input to the system the arrangement 
of Figure 3.13 was used to generate a coloredinput signal. 

The plant used is that of System 1. The test input 1s a 
white Gaussian noise source followed by a single pole filter, 


with the pole position varied from -0.5 to 0.5 for different 
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test runs. The filter was operated under three conditions 
of preprocessing; 
1) No preprocessing. 
2) Inputs U(w) and Y(w) normalized as in (3=75c) 
Soe Multiplication of the input data by a Hamming window 
and then normalizing in the frequency domain as in 
s= 756). 
Figure 3.14 shows two typical plots of the first stage 
lattice parameters. The horizontal axis is the pole position, 


Pere tne input filter of Figure 3.13. The vertical axis is 
(1) 


the value of the lattice parameter, in this case Ky and 
em for cases 1) and 2). For all first stage parameters, 


case 3) results in essentially identical results to case 2). 
It is clear that the preprocessing has made the lattice 
parameters independent of the input filter characteristics. 


Figure 3.15 shows two similar plots for second stage 


~—( 2) ee 
Ki4 and Ky 5 : 


effective, and case 3) is essentially identical to case 2). 


(3) 
ie 


Menat of case 1) and 2). Here the parameter shows some 


parameters Again the preprocessing is 


Figure 3.16 shows parameter k The first plot is 
Bemsitivity to the pole position, though much less than the 
unnormalized input signal. The second plot is for cases 
1) and 3). The windowed and compensated lattice still shows 
no sensitivity to the input filter pole position. Of the 


eight parameters, this one showed the most variation with the 


Pere position. 
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Here the 
parameter calculated without windowing varies over an even 
Wider range than the parameters calculated without any 
preprocessing. Several other of the parameters of the fourth 
stage do the same. The parameters of the windowed and com- 
Bensated jlattice, case 3), are all still insensitive to the 
miput pole position. 

It is apparent from this experiment that input signal 
preprocessing is effective in compensating for the effects 
of a nonwhite input signal, and results in a consistent 
set of lattice parameters for identification. The pre- 
processing should consist of windowing of the input data 
(to eliminate the Gibb's phenomenon) and normalization to 


whiten the signal. 


Meer lLlbR SYNTHESIS WITH THE ARMA FREQUENCY DOMAIN LATTICE 

So far the lattice filters have been developed as models 
for systems with measurable input and output. It is possible 
to define the input signal to the filter as a unit impulse 
so that the lattice models a desired response in the frequency 
or time domain. In this application the recursive in order 
nature of the lattice solution becomes of interest. The 
Ffattice can be expanded in order until a suitable fit of the 


desired specification is obtained. 
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Consider what happens when the frequency domain ARMA 
lattice is used to model a system when the excitation of that 
system is a unit impulse. The configuration is shown in 
Figure 3.18. U(n), the spectrum of the input to the system 
under test, 1s a flat spectrum of unit magnitude, exhibiting 
the previously discussed desirable property of equal energy 
across the spectrum. The output spectrum of the system Y(n) 
is exactly H(n), the discrete Fourier transform of the 


impulse response of the system under test. 
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Figure 3.18 Frequency Domain Modeling 





The lattice parameters obtained as a result of this 
analysis could be used to construct a synthesis model. This 
synthesis model would exhibit the best fit, -in a least mean 
Square sense, that a filter of its order could obtain. The 
order of the lattice can be increased to obtain an arbitrarily 
small MSE. When designing to a specification it is often 
desired that some portions of the frequency characteristics 
of the resultant filter be weighted more heavily than other. 
This too can be accomplished with the ARMA lattice. U(n) and 
Y(n) can both be multiplied by a weighting function which 
provides preemphasis to those portions of the spectrum for 
which a more accurate fit is desired. 

im filter Sumuneene Exanp les 

To illustrate this capability it was postulated that 


a digital filter was desired with a magnitude 
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JH(n)| = 0.495 cos (yap 





n) +.505, O<n<1023 


This desired magnitude was transformed to a minimum phase 
Spectrum using a digital approximation of the Hilbert trans- 
form [14]. This spectrum was used as the input Y(n) to the 
frequency domain lattice, with a unit real spectrum applied 
as the U(n) input. The magnitude of the first and second 


order models of this desired spectrum are shown in Figure 


3.19. The second order model magnitude appears to match 
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the desired magnitude, and as shown by the pole zero plot, 
Pecure 3.19c, 1S Stable and minimum phase. The digital 


meamsrer function of this second order filter is 


il 2 


+0.202512z 
ay, 


0.302504+0.495032z. 
H(z) = =a 
1+#0.00004z ~+0.000005z 


As another example the frequency domain lattice is 
used to find a digital equivalent of the analog low pass 
Gircuit of Figure 3.20 which exhibits 40 dB attenuation at 
5 kHz. This circuit was analysed using the SPICE program, 

a circuit analysis program. This program produced an AC 
analysis of the circuit, showing the complex output when the 
circuit is driven by a unit magnitude signal for frequencies 
from DC to 10 kHz, the equivalent of a Bode plot of the 
@imeuwit transfer function. 

This complex spectrum was applied to the frequency 
domain lattice as in the previous example. The third, 
memeten and fifth order solutions are shown in Figure 3.21. 


Mme transfer function of the fifth order solution is 
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Figure 3.20 Analog Low Pass Filter 





The pole zero and phase plots of the fifth order 
solution are shown in Figure 3.21. There is an additional 
pole not shown on the plot that is located at 3=-4.9. ~The 
Phase plot, Figure 3.2le, shows that the analog circuit has 
a phase at the Nyquist frequency of approximately m/4. Of 
mBemmese the digital filter cannot match this phase, it must 
be either zero or m at the folding frequency. It is 
believed that this discontinuity in the phase is the cause 
of the digital filters non minimum phase characteristic. 
The filter is stable, however, and does meet the desired 


magnitude specification. 
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IV. AUTOREGRESSIVE LATTICE PARAMETERS CALCULATED 


UStiGelheecORREDATION MATRIX 


In the foregoing chapters several methods of solution of 
the AR and ARMA normal equations have been discussed. These 
BPemations can be solved by matrix methods but for large 
order problems this 1s a computationally expensive method. 
If the available data is windowed to ensure a Toeplitz and 
symmetric correlation matrix the equations can be more 
efficiently solved using Levinson's algorithm. From 
Levinson's recursion the lattice structure can be imple- 
mented but requires more calculations because it is necessary 
to update prediction error signals at each stage. Makhoul 
[10] developed an algorithm that implements the single 
channel AR lattice but does not require the generation of 
these error signals, and thus is more efficient than other 
methods. 

In this chapter a new implementation is developed for 
both the single channel and the two channel AR lattice that 
implements the lattice structure without the expense of 
Semeulating error signals. In the final formulation there 
is no requirement that the correlation matrix be Toeplitz 
or symmetric, removing constraints on how data is windowed 
and allowing solution of the problem when the system under 


analysis is not stationary. However, since the formal proof 


val 





depends upon the lattice formulation which assumes symmetric 

Me -vezmadtrices, two SOlutions are presented, one 

assuming Toeplit symmetric correlation matrices, the other 
Smeteing this restriction. This algorithm maintains the 
lattice advantages of maximum entropy and robustness that 

have made the lattice an attractive analysis structure. It 
requires less storage than conventional lattice implementations 
and can be extended efficiently-~to multichannel structures as 


are used for nonlinear systems analysis. 


See LHE SINGLE CHANNEL AUTOREGRESSIVE LATTICE 
In the single channel autoregressive lattice structure 
of Figure 4.1 the forward and backward error signals for 


each stage are updated using the equations from Section II.B. 
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1) 6) and SG) are sums of weighted 


Since the quantities e 
present and past values of the input y(k), the numerator and 
denominator of (4-3) can be expressed as the weighted sum of 
correlations of the input vector een algorithm based upon 
this observation can be expected to be computationally more 
efficient than previous algorithms if the correlations of the 
input signal are known since it is not necessary to update 


the error sequences at each lattice stage. 


From Figure 4.1 the error signal e™ (x) can be written 
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et? is the summation of the gain factors of all transmission 


a 
eeens with delay i from the input to 2D) Oa A recursive 
formulation for we 1s now developed. Equation (4-4) can 
be written 
(n) (ny? (n) 
= (k) = y (k) w (4-6) 
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The weight vector associated with e (k) ore (k) can 


be found recursively. Substituting (4-6) and (4-11), the 
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backward error recursion, results in 
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yielding the third stage weight vectors. 


The numerator of equation (4-3) can now be rewritten 


using the definition of the error signals in terms of the 


weight vectors 
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[ments Seen that the numerator of (4-3). required for calcu- 
Wation of oe can be calculated with no presumption of 
symmetry or requirement that the correlation matrix be 
Toeplitz. The lattice method may still be used if the 
Slgnals under analysis are not stationary or estimates of 
the correlation matrix are available that do not result ina 
Toeplitz and symmetric matrix. However, if the correlation 
matrix of (4-20) is Toeplitz and symmetric the matrix 


multiplication may be written as summations. 
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These summations are recognized as convolutions of the 
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weight vector w , or the correlation of the forward and 
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of two vectors, one with elements of Nee mie eenier of. 02 
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The denominator of equation (4-3) could be found using 
Similar arguments, but generally offers no advantage to 


previously described methods. The recursion 
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(4-20) or (4-26) and (4-29). 
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1. Implementation of the Single Channel Algorithm 


The lattice can be implemented using correlation 


weights by applying (4-16), (4-13), (4-27), (4-20) or 


(4-28), and (4-29) and (4-30). A flowchart of the procedure 


is presented as Figure 4.2. 


For an mth order analysis the procedure will be: 
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Gee. 


correlation vector vr . For nonstationary 


systems or infinite unwindowed data find the 
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Find w For the first stage of analysis the 
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later stages (4-13) gives the recursive solution 


TOO) 


for w 
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vector and either a correlation vector r or 
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the correlation matrix R 
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Pigure 4.2 Flow Chart for mth Order Analysis 
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Cn) Co) 


co) Pamyaien From (4-28), For n=0, P = 8 nye 


otherwise use the recursion of (4-28). 
CS) bat l) me meweroumad from (4-30). 
(6) If higher order analysis is desired, increment 
n and return to (2). 
Figure 4.3 gives a comparison of the number of 
multiplications required to calculate reflection factors for 
a seven stage lattice that uses 256 data points, and assumes 


a stationary system with the data windowed to ensure a 


Toeplitz structure. 


CONDITIONS: 7 stages 
256 data points 
Data windowed for Toeplitz and symmetric 
correlation matrix 
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Eee oat, TWO CHANNEL LATTICE 

The two channel lattice developed in II,C can similarly 
be redefined in terms of correlation weights, the details of 
which are presented in Appendix A. The following is a 
summary of the steps and equations necessary to implement 
the algorithm. 

To implement the two channel algorithm it is necessary 
to calculate eight weighting vectors, two for each of the 
four forward and backward error signals. These can be found 


Weame the initial conditions 
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The vectors of (4-33) and (4-34) each have two elements 
which are also vectors, and the matrix multiplication is 


Peencred as a (2x2) matrix times a (2x1) vector. 


Sa can be found using these weight 
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A flow chart of this implementation is presented as 


Peeure 4.4, 
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Figure 4.4 The m-order Two Channel Implementation 
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iPecorrelattonenatreice Simulations 

To verify the algorithms presented in the forgoing a 
FORTRAN program was written. This program implements the 
flowchart of Figure 4.4. The correlation matrices generated 
by this program are both Toeplitz and symmetric. The 
examples here are the same two systems used in Chapter III.E. 
and represent a run which is typical of the ensemble of tests 
that were made. The distribution over the ensemble was minor 
Since measurement noise was not considered. The fourth order 
model of System 1 is shown with Gaussian elimination using 
256 points, in Figure 4.5. Figure 4.5 shows the results 
using 4096 points. Figure 4.5 shows the same comparison of 
eeyetem 2 using 128 data points, and Figure 4.8 is System 2 
modeled with 1024 points. In almost all cases the results of 
the correlation method are essentially identical to the 
results of the Gaussian elimination method. This should not 
be unexpected. Both matrix manibdulation and correlation 
lattice methods use identical data, the auto and cross 
correlations of the input data signals. Both methods window 
the data identically. Both methods are solving the same 
normal equations, though through differing algorithms. 

The number of multiplications saved by the correla- 
tion lattice, when compared to the conventional lattice 


implementation, 1S approximately 8p, where p is the number 
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of data points. This is because the error signals are not 
multiplied by the lattice parameter matrices ee and mm 


in the correlation implementation. 
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Vc weCONCEUSIONS AND OPEN QUESTIONS 


The objective of this dissertation has been to develop 
efficient methods of implementing the AR lattice. Two 
particular applications for system identification and 
modeling have provided the motivation for this research. 
Mime first is system fault identification. Computationally 
efficient algorithms must be implemented if ie fperuleleate 
analysis is expected to be done on a real time system. 
Adaptive lattice structures require only moderate computation 
for each data point and have often been implemented in real 
time. The block structures, such as the conventional time 
domain or frequency domain lattices, require greater 
mambpers of calculations in bursts for each block of data. 
These structures could be implemented using both parallel 
and pipeline procedures: parallel calculation of the 
correlations of each channel and pipelined calculation of 
the stages of the model. 

The correlation lattice structure also encourages 
parallel processing implementations. In this structure 
calculation of the correlation matrix requires the greatest 
computational effort. These correlation estimates can be 
performed in parallel as the data becomes available. 

The second application considered is the modeling of 


Moniinear systems. For this application the foregoing 


ils) 





issues are of even greater importance. Hybrid signals, 
Peeagucts and cross products of the nonlinear system input 
and output, are generated and the resultant signals applied 
as the input of a multichannel lattice. In this case the 
computational and storage savings of the correlation lattice 
become of major significance and allow the implementation of 
more complex models. It has been found that the two channel 
lattice when applied to ARMA modeling can be simplified 
because of the relationships that exist between the input 
and the output of the system under test. The lattice can 

be simplified further when the system driving signal is white. 

Some of the constraints found are given physical inter- 
pretations that relate to the system under analysis. There 
has been no such interpretation for the relationship of 
the lattice parameters by a factor of the DC gain of the 
system being modeled. It is felt that this simple relationship 
should have some physical interpretation related to the 
Sseem under test. 

The lattice was redefined using the frequency domain 
“pepresentation of the input data. This was found to be 
useful because operations in the frequency domain allowed 
for normalization of the input signals to whiten their 
spectra. This whitened spectra provides a standard input 


that allows lattice parameters to be used for identification 
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and ensures the previously mentioned relationships of the 
lattice parameters. 

There are two normalization procedures proposed. The 
first requires the averaging of ensembles of the inputs to 
obtain an expectation of the energy spectrum of the input. 
The second method uses the spectrum of a single sample 
Spectrum as a normalization function. Either of these 
methods will fail if the energy at any single spectral 
harmonic is zero. 

It 1S not understood how either of these methods affects 
the accuracy of the resultant model. Experiments comparing 
normalization methods, with and without spectral smoothing 
introduced by windowing, have been inconclusive. While some 
examples have shown normalization can improve model accuracy 
over an unnormalized analysis, other examples show the 
opposite. It is postulated that the dynamic range of the 
input spectrum and normalization spectrum are of major 
Slgnificance, but there has been no analysis of the 
procedures involved. 

It is also found that this frequency domain lattice can 
be used to synthesize a digital filter. The desired 
Frequency domain spectrum is provided as the input of the 
lattice, and recursive in order lattice algorithm provides 


an LMS lattice model which meets the desired specification. 


igs! 





Lastly it is found that the lattice can be reformulated 
in terms of correlation weights of the input signals. This 
is efficient because it becomes unnecessary to generate 
error signals within the lattice. This results in both 
computational and storage efficiencies when compared to 
conventional lattice implementations. Though developed for 
the two channel lattice, the procedure can be extended to 
higher dimensioned lattices where the efficiencies will be 
of even greater significance. | 

Inherent in the derivation of Levinson's algorithm is 
mie requirement that the input signals correlation matrix be 
symmetric Toeplitz, requiring that the signals be stationary 
and correlations estimates are from windowed data. When the 
lattice structure is derived from Levinson's algorithm the 
requirement of windowing the data is displaced by the 
lattice structure, which forces upon the data signals a 
different but nevertheless present window function. The 
symmetric Toeplitz nature of the correlation matrix is not 
an issue with the lattice structure. 

When the lattice is reformulated in terms of correlation 
weights, 1t 1S again necessary to generate a correlation 
matrix. It does not appear in the derivation that there is 
a requirement for a symmetric Toeplitz matrix but it is also 
not understood how a nonsymmetric Toeplitz correlation 


matrix will affect the accuracy of the solution. This 
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formulation can be expected to generate errors in the 
solution, and it is not clear how they will propagate 
through the correlation lattice algorithm. This is an issue 


that deserve future attention. 





APPENDIX A 


ine two chemne!l lattice structure of Figure 2.8 can be 
reconfigured for a more efficient implementation using 
correlation weights rather than error signals. The lattice 


equations earlier defined in II.C are repeated here for 


convenience. 
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The result of the multiplication of (A-12) is a 2 element 


Semumn vector. 
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Using this same notation the backward delayed error 
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In the single channel case the backward weighting vector 
was equal to the forward weighting vector inverted in place. 
ents arose because forward and backward reflection factors 
are equal in the single channel lattice. In the two channel 
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Using this notation (A=3) may be rewritten by 


Meecrtucing (A-12) and (A-17), adjusted for proper order, 
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Comparing (A-19) with (A-12) and adjusting for order (ntl), 
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This provides a recursive solution for the forward weight 
meectors. 
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and @*? (x) can be written, with a delay added, from (A-2) as 
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Inspection of Figure 2.8 shows that (A-24) is indeed the 
equation of the first stage error signal. 

An analogous procedure can be allowed to develop the 
recursive solution for the backward weights. Substituting 


_ (A-12) and (A-17) into (A-4) yields 
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Gomparins with (A=17) and adjusting for order (ntl) reveals 
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With the weight vectors for all stages established, the 
reflection factors may be found in terms of these weights. 
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The center matrix is dimensioned (2x2). Its four elements 
are each a submatrix of dimension (2(n+2)x2(n+2)). Each 
submatrix is premultiplied by a 2(n+2) element row vector 
and postmultiplied by a 2(nt2) element column vector. The 
result of these multiplications will be a scalar value for 
each of the four elements of an 


The order of these vector multiplications may be 
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Implementation is discussed in Chapter IV.B. 
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